{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load required packages\n",
    "using Gadfly\n",
    "using Colors\n",
    "using ODEInterface\n",
    "using HDF5\n",
    "@ODEInterface.import_huge\n",
    "loadODESolvers();\n",
    "\n",
    "# Number of subdivisions of the rope\n",
    "global n = 40;\n",
    "\n",
    "# Define the system of ODEs\n",
    "function rope(t,x,dx)\n",
    "    n2 = n*n; # n^2\n",
    "    n3by4 = convert(Int64,3*n/4); # 3*n/4\n",
    "    \n",
    "    # Force in x-direction\n",
    "    Fx = 0.4;\n",
    "    # Force in y-direction\n",
    "    Fy = cosh(4*t-2.5)^(-4);\n",
    "    \n",
    "    # Compute required matrices\n",
    "    c = -cos(x[1:n-1]-x[2:n]);\n",
    "    cDiag = [1;2*ones(n-2);3];\n",
    "    C = spdiagm((c,cDiag,c),(-1,0,1));\n",
    "    \n",
    "    d = -sin(x[1:n-1]-x[2:n]);\n",
    "    D = spdiagm((-d,d),(-1,1));\n",
    "    \n",
    "    # Compute the inhomogeneous term\n",
    "    v = -(n2+n/2-n*[1:n;]).*sin(x[1:n])-n2*sin(x[1:n])*Fx;\n",
    "    v[1:n3by4] = v[1:n3by4] + n2*cos(x[1:n3by4])*Fy; \n",
    "    \n",
    "    w = D*v+x[n+1:2*n].^2;\n",
    "    u = C\\w;\n",
    "    \n",
    "    # Write down the system\n",
    "    dx[1:n] = x[n+1:2*n];\n",
    "    dx[n+1:2*n] = C*v + D*u;\n",
    "    \n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Initial Conditions\n",
    "t0 = 0.0; T = 3.723; x0=zeros(2*n);\n",
    "\n",
    "# Compute the \"reference solution\"\n",
    "opt = OptionsODE(OPT_EPS => 1.11e-16,OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "OPT_RTOL => 1e-16,OPT_ATOL=>1e-16);\n",
    "(t,x_ref,retcode,stats) = odecall(dop853,rope,[t0, T], x0, opt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = zeros(length(t),n+1);\n",
    "y = zeros(length(t),n+1);\n",
    "l=1/n;\n",
    "\n",
    "for j = 1:length(t)\n",
    "    for i=1:n\n",
    "        x[j,i+1] = x[j,i]+l*sin(x_ref[j,i]);\n",
    "        y[j,i+1] = y[j,i]+l*cos(x_ref[j,i]);\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h5write(\"ropePlotx.h5\",\"x\",x)\n",
    "h5write(\"ropePloty.h5\",\"y\",y)\n",
    "h5write(\"ropePlott.h5\",\"t\",t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
